Contents

Model vs. satellite sea ice thicknessΒΆ

In this notebook, we compare remote sensing-derived sea ice thickness from ICESat-2 and model-derived sea ice thickness from the Pan-Arctic Ice Ocean Modeling and Assimilation System (PIOMAS, Zhang and Rothrock, 2003)

import xarray as xr # For working with gridded climate data 
import pandas as pd
import numpy as np
import itertools
from utils.read_data_utils import read_book_data # Helper function for reading the data from the bucket
from utils.misc_utils import restrictRegionally # Region restriction
from utils.plotting_utils import compute_gridcell_winter_means, interactiveArcticMaps, interactive_winter_mean_maps, interactive_winter_comparison_lineplot # Plotting

# Plotting dependencies
import cartopy.crs as ccrs
from textwrap import wrap
import hvplot.xarray
import holoviews as hv
import matplotlib.pyplot as plt
from matplotlib.axes import Axes
from cartopy.mpl.geoaxes import GeoAxes
GeoAxes._pcolormesh_patched = Axes.pcolormesh # Helps avoid some weird issues with the polar projection 
%config InlineBackend.figure_format = 'retina'
import matplotlib as mpl
mpl.rcParams['figure.dpi'] = 150 # Sets figure size in the notebook

# Remove warnings to improve display
import warnings 
#warnings.filterwarnings('ignore')

1) Read in dataΒΆ

Below, we’ll read in the data and restrict it to the Inner Arctic region. For the PIOMAS sea ice thickness, we’ll mask out the pole hole and set 0 thickness to nan to match the ICESat2 data; the PIOMAS data sets anywhere with no sea ice to nan, whereas PIOMAS sets anywhere with no sea ice to 0.

book_ds = read_book_data() # Read/download the data 
book_ds = book_ds[["ice_thickness_smoothed","piomas_ice_thickness"]] # Grab data variables of interest
book_ds = restrictRegionally(book_ds, regionKeyList=[10,11,12,13,15]) # Restrict data to the Inner Arctic

book_ds["piomas_ice_thickness"] = book_ds["piomas_ice_thickness"].where(book_ds["piomas_ice_thickness"]!=0, np.nan) # Set 0 --> nan 
book_ds = book_ds.where(book_ds.latitude<=88, np.nan) # Mask out pole hole

years = [2018,2019,2020] # Years over which to perform analysis

Regions selected: Inner Arctic

2) Map mean winter differencesΒΆ

Here, we’ll compute winter means for PIOMAS and ICESat-2 sea ice thickness. Then, we’ll use the interactiveArcticMaps function to generate maps of the winter means, along with a comparison map showing the gridcell differences.

is2_winter_means = compute_gridcell_winter_means(book_ds["ice_thickness_smoothed"], years=years, start_month="Nov")
pio_winter_means = compute_gridcell_winter_means(book_ds["piomas_ice_thickness"], years=years, start_month="Nov")

for i in range(len(years)): 
    data1 = is2_winter_means.isel(time=i)
    data2 = pio_winter_means.isel(time=i)
    pl1 = interactiveArcticMaps(data1, vmin=0, vmax=4, colorbar=False, title=data1.time.values.item())
    pl2 = interactiveArcticMaps(data2, vmin=0, vmax=4, colorbar=True, clabel="sea ice thickness (m)", title=data2.time.values.item()) 
    diff_pl = interactiveArcticMaps((data2-data1), vmin=-1.5, vmax=1.5, cmap="coolwarm", clabel="gridcell difference", title="Gridcell difference") 
    data_pl = (pl1+pl2+diff_pl).opts(hv.opts.Layout(shared_axes=False, merge_tools=True))
    data_pl.opts(title="Winter "+str(years[i])+"-"+str(years[i]+1))
    display(data_pl)
    print("\n")

3) Mean monthly sea ice thickness comparisonΒΆ

Here, we’ll generate an interactive lineplot of the mean monthly sea ice thickness for PIOMAS and ICESat-2.

# Compute monthly means 
mean_monthly_is2 = book_ds["ice_thickness_smoothed"].mean(dim=("x","y"), keep_attrs=True)
mean_monthly_pio = book_ds["piomas_ice_thickness"].mean(dim=("x","y"), keep_attrs=True)

# Generate individual plots
is2_lineplot = mean_monthly_is2.hvplot.line(grid=True, label="ICESat-2",frame_width=600, frame_height=300) * mean_monthly_is2.hvplot.scatter(marker='s') # Overlay scatter plot to add markers
pio_lineplot = mean_monthly_pio.hvplot.line(grid=True, label="PIOMAS") * mean_monthly_pio.hvplot.scatter(marker='o') # Overlay scatter plot to add markers

# Combine plots and display
comb_plot = (is2_lineplot*pio_lineplot).opts(hv.opts.Layout(shared_axes=True, merge_tools=True))
comb_plot.opts(title="PIOMAS and ICESat-2 sea ice thickness (Sep 2019 - Apr 2021)")
comb_plot

A quick comparison between PIOMAS and ICESat-2 ice thickness shows that PIOMAS sea ice thickness is consistently higher than ICESat-2, except for in the beginning of each winter, where ICESat-2 is higher. Below, we’ll use the interactiveArcticMaps function to visualize any interesting patterns in the differences for September 2019, the month with the large difference between PIOMAS and ICESat-2 ice thickness (0.9 meters vs. 1.38 meters).

# Grab data for a given month
date = "Sep 2019"
pio = book_ds["piomas_ice_thickness"].sel(time=date)[0]
is2 = book_ds["ice_thickness_smoothed"].sel(time=date)[0]

# Generate maps and display combined figure
pl_pio = interactiveArcticMaps(pio, vmin=0, vmax=4, colorbar=False, title="PIOMAS ice thickness ("+date+")") 
pl_is2 = interactiveArcticMaps(is2, vmin=0, vmax=4, title="ICESat-2 ice thickness ("+date+")") 
diff_pl = interactiveArcticMaps((pio - is2), vmin=-1.5, vmax=1.5, clabel="gridcell difference", cmap="coolwarm", title="Gridcell difference") 
data_pl = (pl_pio+pl_is2+diff_pl).opts(hv.opts.Layout(shared_axes=False, merge_tools=True))
data_pl